Solutions

Questiohs

1. define projection.
2. do some editing.
3. intersect with an existing shapefile.
4. calculate area of obtained polygons

Packages

library(tidyverse)
library(sf)
library(tmap)
library(leaflet)

Data

ucsfi <- read_sf("data/ucsfi.shp")  # %>% `st_crs<-`(3857)  # 3857 is Web mercator
brazil0 <- read_sf("data/BRA_adm/BRA_adm0.shp")

Projection

In the above code chunk projection is commented out because tmap appears to work beset without the projection assigned. We’ll do that later. assigned with st_crs() using a somewhat non-standard syntax. However, It seemed to be the only way I could add a projection to this particular shapefile. Below are some other syntax examples that have worked for me with other shapefiles. Of course there are various ways to add projection and since this is R, there’s more than one way to accomplish the goal.

sf::st_crs(3857)
sf::st_transform(3857)
ggplot2::coord_sf(crs = 3857)   # must use the development version of ggplot2

Notes on Projection

From Brazillian CRS

• 4989
• 4674
• 5527
• 4618

Check Projection

st_crs(ucsfi)

Coordinate Reference System: NA

print("ucsfi above ^^^  --------- brazil0 below")

[1] "ucsfi above ^^^  --------- brazil0 below"

st_crs(brazil0)

Coordinate Reference System:
  EPSG: 4326 
  proj4string: "+proj=longlat +datum=WGS84 +no_defs"

# or check the projection for an individual polygon ...
# ucsfi$geometry %>% head(1) 

First View With Projection

Data Preview

as_tibble(ucsfi)

Simple feature collection with 147 features and 13 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -73.99068 ymin: -32.92342 xmax: -32.35011 ymax: 5.272257
epsg (SRID):    NA
proj4string:    NA

Wrangle Data

I want to use a numeric data type as my “fill”

ucsfi <- ucsfi %>% 
  mutate(GID7 = as.numeric(GID7))

First Look with tmap

tmap_mode("plot")

tmap mode set to plotting

tm_shape(brazil0, projection = 3857) + tm_fill() + tm_text("NAME_0", auto.placement = FALSE) +
  tm_shape(ucsfi, projection = 3857) + tm_fill("GID7", palette = rev(sf.colors()))

Currect projection of shape ucsfi unknown. Long-lat (WGS84) is assumed.

ggplot

ucsfi %>% 
  mutate(ANO_CRIA6 = as.numeric(ANO_CRIA6)) %>% 
  ggplot() +
  geom_sf(aes(fill = ANO_CRIA6, color = ANO_CRIA6)) 

NAs introduced by coercion

Look 3 – Leaflet

MapPalette <- colorQuantile(palette = "viridis", domain = ucsfi$GID7, n = 7, reverse = TRUE)

ucsfi %>% 
  leaflet(width = "100%") %>% 
  addProviderTiles( provider = "CartoDB.Positron") %>% 
  addPolygons(popup = ~ ID_UC0,
              stroke = FALSE,
              smoothFactor = 0,
              fillOpacity = 0.7,
              color = ~ MapPalette(GID7)) 

#+
#  addLegend("bottomright", 
#              pal = MapPalette, 
#              values = ~ GID7,
#              title = "GID7",
#              opacity = 1)